1. Field of the Invention
The present invention comprises methods for generating surface-consistent residual statics. In particular, the invention relates to methods for calculating static corrections for seismic data by minimizing an objective function using global optimization thereby to maximize a corresponding stack power. The stack power is decoupled from the overlapping between different seismic traces into several one-dimensional problems. The surface-consistent residual statics that minimize the objective function are determined by using a Stochastic Pijavskij Tunneling technique to eliminate regions where the global minimum of the objective function is unlikely to exist so that the global minimum of the objective function can be reached without determining the objective function in these eliminated regions.
2. Background Art
The Problem
Seismic exploration is a predominant geophysical activity conducted to find commercial accumulations of oil, gas, or other minerals, to study the nature of the Earth for the foundations of roads, buildings, dams, tunnels, nuclear power plants, and other structures, and to search for geothermal areas, water resources, archeological ruins, etc.
In seismic exploration, seismic waves or signals are generated by shots from one of several types of energy sources and detected by arrays of sensitive devices or receivers called geophones or hydrophones. The most common measurement made is of the travel times of seismic waves, although attention has been increasingly directed to the amplitude of seismic waves or changes in their frequency content or wave shape.
The measurement of seismic signals may be seriously interfered by external factors such as wind, vehicular and pedestrian traffic, microseisms and at sea, ship noise and environmental noise due to marine life. Moreover, solid friction in the Earth as well as seismic waves scattering attenuate the measured signals. Thus, a plurality of traces from different shots resulting from the interaction of a plurality of cooperating sources and receivers, but having the same subsurface incident point, are gathered and combined together by use of various algorithms well known to those skilled in the geophysical art. Multiple coverage tends to destructively attenuate random noise but enhance coherent reflection data. One such method is referred to as a Common Mid Point ("CMP") stacking which will now be illustrated with reference to FIG. 1.
In FIG. 1, sound waves generated by acoustic sources S.sub.1 -S.sub.3 located on surface respectively are detected by receivers R.sub.1 -R.sub.3, after sound wave reflection from interface 14 along ray paths 20 and 30, 22 and 28, 24 and 26. The number of the energy sources can be different from that of the receivers. The ray paths shown all converge at a common mid point 16. The surface separation between the respective sources and correspondingly numbered receivers is termed the offset. Correction for angularity (also called Normal Moveout or "NMO") to the respective ray paths by stacking them together would collapse them to zero-offset as represented on the surface by S.sub.0 /R.sub.0 to form two-way travel path 18 impinging on the common mid point 16. In this case of zero dip, mid point 16 also forms a common depth point as well. Because all of the ray paths are incident on the same subsurface point, here CMP 16, the traces can be stacked validly to enhance the signal-to-noise ratio.
Near-surface lateral velocity variations and surface elevation changes create travel-time variations that may be approximated by surface-consistent static time shifts. In FIG. 1 there are shown two subsurface earth layers I and II, separated by line 12, characterized respectively by acoustic propagation velocities of V.sub.I =1500 meters per second (m/s) and V.sub.II =3100 m/s. P.sub.si and P.sub.rj identify those portions of the near-vertical ray paths that traverse the variable-thickness upper low velocity layer I, after ray bending due to Snell's law of refraction at the interface 12, beneath sources S.sub.i and receivers R.sub.j. Because of a longer combined path length of P.sub.s3 +P.sub.r3 through the low velocity layer I, the total travel time along ray path 20, 30 from S.sub.3 to R.sub.3, after application of NMO, will still be longer than the total travel time along either ray path 22, 28 from S.sub.2 to R.sub.2 or ray path 24, 26 from S.sub.1 to R.sub.1, corrected for NMO, where the combined ray path P.sub.s2 +P.sub.r2 or P.sub.s1 +P.sub.r1 through the low velocity layer I is shorter. Note that FIG. 1 gives an example where the velocity in the low velocity Layer I has a constant volume. In reality, this low velocity Layer I, which is usually 4 to 50 meter thick, is characterized by seismic velocities that are not only low (usually between 200 and 2,000 m/s), but at times highly variable. Proper statistical methods need to be used to estimate the low velocity layer. It is known in the art, therefore, that even though the low-velocity segments of the total ray path are relatively short, travel-time differences are not negligible due to the very low velocity at or near the surface I. The respective time delays s.sub.i and r.sub.j due to a variable-thickness low velocity layer, here the surface layer I, are defined as the surface-consistent statics which must be applied as corrections to the reflection travel times prior to stacking for maximizing inter-trace reflection coherency. The statics are termed surface-consistent because they are due to irregularities of the near-surface low velocity layer.
In the example of FIG. 1, only a few of sources and receivers are shown for simplicity. The relative static corrections can be easily obtained for this simple, constant velocity case. In reality, calculating surface-consistent residual statics correction to compensate for time shifts in stacked many thousands of seismic traces, many of indifferent or poor quality, recorded from interfacing between up to thirty cooperating seismic sources and receivers poses a challenge to the art.
The Prior Art
Since the earliest days of seismic exploration, geophysicists have recognized the need to correct for the low velocity in the weathered and unconsolidated sediments near the earth's surface. The data processing procedure has been described by Yilmaz, in a book entitled "Seismic Data Processing," published by Society of Exploration Geophysics, Tulsa, Okla., 1987, Marsden, in Static Corrections--a Review, The Leading Edge, 12, Part 1, pp. 43-49, Part 2, pp. 115-120, and Part 3, pp. 210-216, 1993, and Sheriff and Geldart, in Exploration Seismology, Second Edition, Cambridge University Press, Cambridge, England, 1995. The first corrections for elevation and low velocity are field statics. A reference level is determined that is below the low velocity layer (LVL) and field statics move the sources and receivers to the reference level. CMP gathers are used to generate a set of preliminary velocity picks that are used to calculate NMO corrections. Residual statics corrections are calculated using the corrected data. The process is repeated until the results converge.
The conventional method for calculating residual statics corrections was developed by Taner, et al., in "Estimation and correction of near-surface time anomalies," Geophysics, 39, pp. 441-463, 1974 and Wiggins, et al., in "Residual statics analysis as a general linear inverse problem," Geophysics, 41, pp. 922-938, 1976. The time delays caused by the passage of seismic signals through the LVL should depend on path. After the NMO corrections, it is assumed that all of the paths are vertical and estimate a single time delay that is "surface consistent" (each source and receiver location has a time delay that does not depend on the wave path). The time delay, denoted by t.sub.srk, for a trace that follows a path from a source (s) to a receiver (r) via a common midpoint (k) is determined by maximizing the cross correlation between the trace and the CMP gather. The total time delay has four components: the source static (S.sub.s), the receiver static (R.sub.r), the two way travel time from the reference level to a reference subsurface reflector (.GAMMA..sub.k), and a residual NMO correction [M.sub.k (X.sub.sr.sup.2)] where X.sub.sr is the distance from the source to the receiver, as shown in equation (1): EQU t.sub.srk =S.sub.s +R.sub.r +.GAMMA..sub.k +M.sub.k (X.sub.sr.sup.2). (1)
Although the total time delay (t.sub.srk) is an independent value for each trace, the four types of parameters in the total time delay equation each occur in many traces. Thus, the parameters are overdetermined (i.e., there are many more equations than parameters). For example, for a problem with 4776 traces, 100 shots, 216 receivers, and 423 common midpoints, there are 4776 equations and 1162 parameters. Since the equations are overdetermined, they can only be solved by minimizing the errors (by using the method of least squares).
Taner et al. discovered that the problem is both overdetermined and underdetermined. Although there are many more equations than unknowns, he found sets of nonzero parameters that satisfy the total time delay equation when all of the total time delays are zero. These solutions form a null space and we can add an arbitrary linear combination of the vectors in the null space to any solution of the total time delay equation and still have a solution. Taner et al. added an extra set of terms to the least squares objective function to eliminate the null space solutions.
Wiggins et al. clarified the work of Taner et al. by applying singular value decomposition to find the least squares solution of the total time delay equation, by displaying the eigenvectors for a small problem, and by introducing the term "null space."
Residual statics corrections can be calculated in two steps: use crosscorrelation to estimate the total time delay (t.sub.srk) for each trace and use least squares to find the parameters in the total time delay equation. Ronen et al., in "Surface-consistent residual statics estimation by stack-power maximization," Geophysics, 50, pp. 2759-2767, 1985, proposed a one step alternative approach: stack power maximization. Ronen et al. defined an objective function that measures the correlation between all of the traces in each CMP gather. Changes in the parameters in the total time delay equation cause a time shift for each trace and change the correlation between traces. Ronen et al. searched for parameter values that maximize the stack power.
If synthetic data is created by time shifting a trace, correlation of the two traces will identify the time shift (the correlation will be 1.0 at the proper time shift). For real data, there may be many local maxima in the correlation function. The failure to find the largest local maxima results in a "cycle skip." The stack power function depends on thousands of traces and hundreds of parameters and can have a very large number of local maxima (if there are N parameters and M local maxima in each dimension of the N dimensional space, the logarithm of the total number of local maxima is Nx log M). Most optimization methods will find a local maximum but not a global maximum. A problem with many local maxima requires a global optimization method.
Rothman, in "Nonlinear inversion, statistical mechanics, and residual statics corrections," Geophysics, 50, pp. 2784-2796, 1985, recognized that the residual statics problem was a global optimization problem and proposed to solve the problem using the simulated annealing method. He defined the stack power optimization using two of the four types of parameters in the total time delay equation; his parameters are the static corrections for the sources and receivers. Since the total stack power is a sum of the stack power for each CMP, the CMP term (.GAMMA..sub.k) shifts all of the traces in the stack and does not change the stack power. He argued that the residual NMO correction [M.sub.k (X.sub.sr.sup.2)] was cumbersome to estimate and had a small impact on the stack power. Subsequently, Rothman, in "Automatic estimation of large residual statics corrections," Geophysics, 51, pp. 332-346, 1986, improved his method and applied it to some field data that required large static corrections (up to 200 ms).
DuBose, in "Practical steps toward realizing the potential of Monte Carlo automatic statics," Geophysics, 58, pp. 399-407, 1993, proposed several innovations to improve Rothman's method. The simulated annealing method depends on a parameter: the "temperature." DuBose proposed an algorithm for finding the proper temperature. To apply his method, he measures the changes in the statics from one iteration to the next. If the changes are too small, he raises the temperature. If the changes are too large, he reduces the temperature. Since he performs calculations in the frequency domain, each trace can be stored compactly and time shifting is performed by an efficient multiplication by a complex number. In the time domain, each trace is sampled at a regular period and cross correlation for a time shift that is a fraction of a sampling step requires interpolation. No interpolation is required in the frequency domain.
DuBose assumes that receiver statics are constrained with position along the line such that the change from one receiver to the next is bounded by a value substantially smaller than the largest possible static. He also assumes that any shot that is very close to a receiver station should have its static constrained within another modest limit to be near the static of its receiver station. These assumptions limit the regions in parameter space that need to be searched. He modifies the objective function by performing a running average over a few adjacent stacked traces. To eliminate the null space from the problem, he adds a penalty term to the objective function.
U.S. Pat. No. 5,508,914 issued to Lee disclosed a method for calculating static corrections for seismic data by minimizing an objective function. The surface consistent source and receiver statics that maximize the coherence of reflections are determined by simulating the crystallization process of a system of fictitious classical particles whose interaction is given by the objective function.
Despite many efforts and progresses made in the field, the methods cited above are quite expensive of computing time. None seems to answer the following two inquiries very well: (1) how to find the correct maximum if there are many local maxima, and (2) how to get to the correct maximum with the fewest calculations. Therefore there is a need to develop a more efficient method to estimate the residual statics for massive numbers of seismic traces.
Recently, an improved algorithm for solving global optimization problems was developed at the Center for Engineering Science Advanced Research (CESAR) at the Oak Ridge National Laboratory (ORNL), as disclosed by Barhen, et al., in "TRUST: A deterministic algorithm for global optimization," Science, 276, pp. 1094-1097, 1997. The algorithm is called TRUST (terminal repeller unconstrained subenergy tunneling). The results reported in Barhen, et al. were fruits of a collaborative project to apply TRUST to the residual statics problem using an objective function defined by DuBose, although the objective function was not identical to the function used in DuBose (1993). While the result was promising, it still took a substantial computing time to find a global local minimum. Indeed, it did not always find the globally optimum values for large residual statics problems with 100 or more parameters. This, combined with expensive computing time associated with the current available methods, has prompted further investigation of improvements in calculating the static corrections for seismic data.